% simule.m 

% Simulate the model by feeding generated shocks
% (under full information or incomplete information) 

% then compute and report various results/moments about those simulated data, 


nM = rows(M);
nPI = rows(PI);
nstates = 4;
nk = 5;

ic = 1;   % position of consumption in IR
iy = 7;   % position of output in IR
ii = 9;   % position of investment in IR
iva = 12; % position of value added IR
ispot = 13;
ifor1 = 15;
ifor2 = 16;
ifor3 = 17;
ifor4 = 18;
iinv = 25;
is = 25;  % position of oil inventories in IR



nsample = 250;    % # of periods in each sample
cutoff = 50;      % # of periods cut off from each sample (to reduce
                  % the importance of initial periods)
nreplic = 500;    % # of samples simulated

rng(10);  % so the same stream of "random numbers" is realized every time

myhp = 1;     % 0 if you DON'T want all simulated data to be hp-filtered
              % 1 if you want all simulated data to be hp-filtered
               
kalresults = zeros(nreplic,6); % results from the Kalmn filter : 
                             % position 1 : estimate of phi_tau
                             % position 2 : estimate of sig_p
                             % position 3 : estimate of sig_t
                             % position 4 : frac of var(diff(spot))
                             % explained by permanent shocks according to
                             % Kalman filter estimates
                             % position 5 : value of log l at optimum
                             % position 6 : exit flag from fminunc

                             
% various moments of simulated data
cor1 = zeros(nreplic,nPI);   % correlation between spot price of oil and all other economic variables 
cor2 = zeros(nreplic,nPI);   % correlation of inventories with all other economic variables
sd = zeros(nreplic,nPI);     % standard deviations of all model variables, relative to std(Y)
sdlevel = zeros(nreplic,nPI);     % standard deviations of all model variables, in levels
sdspot = zeros(nreplic,nPI); % standard deviations of all model variables, relative to std(spot)

%matrix of autocorrelation of (HP-filtered) GDP
autocor_y = zeros(nreplic,1);

%matrix of autocorrelation of (HP-filtered) spot price of oil
autocor_spot = zeros(nreplic,1);

%matrix of autocorrelation of (HP-filtered) futures price
autocor_f1 = zeros(nreplic,1);
autocor_f4 = zeros(nreplic,1);
 

 
   
 for j = 1:nreplic

   
   if rem(j,50) == 0
      disp([' Simulation No. ' num2str(j) ])
   end
   
   % create a 'life' of <nsample> periods
   
   IR = zeros(nsample,nPI);
   
   S = zeros(nM,1);    % initial state variables 
                       % equivalent to M*s_t-1 in notes
   shock = zeros(nM,1); %innovations to state variables
   
   z = 0;
   
   stry = zeros(1,1);
   mytry = zeros(nsample,1);
   
   for i = 1:cutoff+nsample; 
         
      temp = randn(nstates,1).*sqrt(diag(QQ));  
      
      shock = [ zeros(nk,1)
                temp 
                zeros(nstates,1) ]; 
      
     
      % augment the state vector with these shocks
      S = S+shock;
      % S is now equivalent to M*s_t-1 + eps_t, in KW terminology
         
      %generate and store corresponding 'flow' variables if <i> is bigger than <cutoff>
      if i > cutoff;
         IR(i-cutoff,:)=(PI*S)';
      end
      
      
      %generate the next "expected" S: 
      S=M*S;
      %Now S stands for M*s_t     
   end
   
   
   

   %--------------------------------------%
   % we now have a sample of simulated data saved in <IR>
   % recall that these are expressed in % deviations from s.s
   
   
   % - Compute the one variable kalman filter decomp as per the empirical
   % part of the paper
   %--------------------------------------%
   
   
   % express spot rates in levels rather than deviations from s.s. 
   spot  =  100*(1+IR(:,ispot)).*po;
   spot;
   
   % compute change 
   data = diff(spot);
   data';
   
   
   % get estimates for phi, sig_tau and sig_P
   
   if(1) % via the Kalmn filter estimates
   
   optz = optimset('fminunc');
   optz = optimset(optz,'Display','off'); 
   optz = optimset(optz,'Largescale','off'); 
   sv1 = [0   0.5    0.0798    0.1556];
 
   % send to Kalman filter to get estimates for phi sig_tau and sigP
   rvll =   @(params)llfpermtempdecompdiff(data',params);
   [rvbout,fval,exitflag] = fminunc(rvll ,sv1(2:end),optz);
   
   % keep results
   sigp = rvbout(1,2)^2;
   sigt = rvbout(1,3)^2;
   phi  = rvbout(1,1);
   vartempcomp = sigt*((1-phi)^2)/(1-phi^2);
   temp  =  sigp/(sigp+sigt+vartempcomp);

   
   %store results from this run
   kalresults(j,1) = phi;
   kalresults(j,2) = rvbout(1,2);
   kalresults(j,3) = rvbout(1,3);
   kalresults(j,4) = temp;
   kalresults(j,5) = fval;
   kalresults(j,6) = exitflag;
   end   
   
   
   % - Compute some relevant moments of these simulated data
   %--------------------------------------%
   
   if myhp == 0  % NO hp filter
       
   cyclic = IR;
   
   else % DO HP filter simulated data
       
   [cyclic,yt] = hpfilt(IR(1:nsample,:),1600);
   
   end
   
   temp = corrcoef(cyclic);
   
   % correlation between spot price of oil and other variables
   cor1(j,:) = temp(ispot,:);
   
   % correlation between inventories and other variables
   cor2(j,:) = temp(iinv,:);
   
   % Volatilities
   temp = real(cov(cyclic));
  
   % compute vector of standard deviations:
   sdtemp = 100*sqrt(diag(temp))';
   
   % save s.d. relative to sd(Y)
   sd(j,:) = sdtemp(1,:)/sdtemp(1,iy); 
   sdlevel(j,:) = sdtemp(1,:);
   
   % save s.d. relative to sd(spot)
   sdspot(j,:) = sdtemp(1,:)/sdtemp(1,ispot); 
   
 
   % first autocorrelation of Y
   temp = xcorr(cyclic(:,iy),1,'normalized');
   autocor_y(j,1) = temp(3);
   
   % first autocorrelation of spot price of oil
   temp = xcorr(cyclic(:,ispot),1,'normalized');
   autocor_spot(j,1) = temp(3);
   
   % first autocorrelation of (one-quarter-ahead) futures
   temp = xcorr(cyclic(:,ifor1),1,'normalized');
   autocor_f1(j,1) = temp(3);
   
    % first autocorrelation of (one-year-ahead) futures
   temp = xcorr(cyclic(:,ifor4),1,'normalized');
   autocor_f4(j,1) = temp(3);
  
 end
 
 
 
 if(1)
 % Report results about various simulated moments
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
disp('               ')
disp([ ' Stylized Facts for simulated data (' num2str(nreplic) ' replications of ' num2str(nsample) ' periods each)' ]) 
disp('*------------------------------------------------------------------------------*')
disp('               ')
   
mvar = mean(sdlevel);
choose1 = [iy ic ii is ispot ifor1 ifor4];
mystring1 = ['GDP                     '
             'Consumption             '
             'Investment              '
             'Oil Inventories         '
             'Spot price of oil       '
             'One-quarter-ahead Future' 
             'One-year-ahead Future   ']; 

table   = [ mvar(choose1)'  ];

disp('   ')
disp('             Table 1: Volatilities ')
disp('*************************************************************************')
disp(['variable (X)      ',  '  std.(X)         ']);
disp('*************************************************************************')
disp([mystring1 num2str(table,'%11.3f')])
disp('  ')




mvar1 = mean(sd);
choose1 = [iy ic ii is ispot ifor1 ifor4];
mystring1 = ['GDP                     '
             'Consumption             '
             'Investment              '
             'Oil Inventories         '
             'Spot price of oil       '
             'One-quarter-ahead Future' 
             'One-year-ahead Future   ']; 

table1  = [ mvar1(choose1)'  ];

disp('   ')
disp('             Table 2: Relative Volatilities (w.r.t. Y) ')
disp('*************************************************************************')
disp(['variable (X)      ',  '  std.(X)/std.(Y) ']);
disp('*************************************************************************')
disp([mystring1 num2str(table1,'%11.3f')])

mvar1 = mean(sdspot);
table1  = [ mvar1(choose1)'  ];

disp('   ')
disp('             Table 2: Relative Volatilities (w.r.t. spot)')
disp('*************************************************************************')
disp(['variable (X)      ',  '  std.(X)/std.(spot) ']);
disp('*************************************************************************')
disp([mystring1 num2str(table1,'%11.3f')])





mvar2 = mean(cor1);
choose2 = [iy ic ii iinv ispot ifor4];
mystring2 = [ ' GDP                    '
              ' Consumption            '
              ' Investment             '
              ' Oil Inventories        '
              ' Spot price of oil      '
              ' One-year-ahead Future  '];
   
table2  = [ mvar2(choose2)'  ];


disp('   ')
disp('             Table 3: Correlations of spot price with other variables ')
disp('*************************************************************************')
%disp(['variable (X)      ','std.(X)  (std) ', 'std.(X)/std.(Y) ',' autocorrelation(X) (std)']);
disp(['variable (X)      ',  '  corr (spot, X) ']);
disp('*************************************************************************')
disp([mystring2 num2str(table2,'%11.3f')])


mvar3 = mean(autocor_y);
mvar4 = mean(autocor_spot);
mvar5 = mean(autocor_f1);
mvar6 = mean(autocor_f4);

mytemp = [table
    mvar4
    mvar5
    mvar6];

disp('   ')
disp('             Table 4: Some Autocorrelations                             ')
disp('*************************************************************************')
%disp(['variable (X)      ','std.(X)  (std) ', 'std.(X)/std.(Y) ',' autocorrelation(X) (std)']);
disp(['variable (X)      ',  '  corr (spot, X) ']);
disp('*************************************************************************')
disp([' GDP                       ' num2str(mvar3,'%11.3f') ])
disp([' Spot price of oil         ' num2str(mvar4,'%11.3f') ])
disp([' One-Quarter-ahead Futures ' num2str(mvar5,'%11.3f') ])
disp([' One-Year-ahead Futures    ' num2str(mvar6,'%11.3f') ])


 end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


 mycount = 0;
 for k = 1:rows(kalresults);
     
     if kalresults(k,6) == 1
         mycount = mycount + 1;
     end
 end
 
  % keep results arrived at with normal convergence
 myoutput = kalresults(kalresults(:,6)== 1,:);
 
 
 % print details if needed
 if(0)
 disp('   ')
 disp('*************************************************************************')
 disp([ ' Number of runs where fminunc converged normally :  '   num2str(mycount)])
 disp('*************************************************************************')
 
 
 disp('   ')
 disp('*************************************************************************')
 disp( ' Results from the Kalman Filter Estimation')
 disp('*************************************************************************')
 disp( '    run       phi_tau  sig_p      sig_tau     frac_p log-lik     exitflag ')
 run = 1:rows(myoutput);
 disp([ run' myoutput])
 end
 
 
 